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A  TIME-DEPENDENT  THREE-DIMENSIONAL  SIMULATION  OF  THE 
EARTH’S  MAGNETOSPHERE:  RECONNECTION  EVENTS 


Introduction 

Tirae  dependent  global  simulations  of  the  earth's  magnetosphere  have 
been  performed  for  several  years.  The  purpose  of  such  a  simulation  has 
been  and  still  is  to  obtain  information  concerning  the  general 
topological  structure  of  the  earth's  magnetosphere  and  to  address 
possible  cause  and  effect  relationships  leading  to  dynamic  phenomena  such 
as  magnetic  substorms.  These  large  scale  calculations  can  also  be  used 
as  boundary  condition  calculations  for  investigation  into  the  relevant 
microphysics  that  produce  small  scale,  but  important  effects.  Finally  it 
is  hoped  that,  as  an  additional  benefit,  large  scale  time  dependent 

simulations  might  provide  further  insight  into  magnetospheric  dynamics, 
so  that  experimental  data  might  be  correlated  with  the  calculation,  and  a 
better  understanding  of  the  earth's  magnetosphere  could  be  reached. 

The  global  simulations  to  be  discussed  below  all  have  one  thing  in 
common.  They  are  time  dependent  solutions  of  the  MHD  equations.  Other 
than  boundary  conditions,  none  of  these  calculations  assumes  a  particular 
structure  for  the  magnetosphere  nor  do  they  have  any  adjustable 

parameters  to  affect  the  structure.  The  number  of  ad  hoc  assumptions  is 
held  to  those  required  for  the  validity  of  the  MHD  equations. 

Customarily  in  a  model,  macroscopic  structure  may  be  manipulated,  e.g. , 
the  magnetopause  boundary  might  be  positioned  at  will.  In  these 
simulations,  only  the  microscopic  structure  -  the  underlying  physics  -  is 
determined.  The  macroscopic  structure  then  either  fits  what  it  is  being 
simulated,  or  it  does  not.  While  questions  of  accuracy  of  solution  are 
always  important,  what  the  simulations  test  is  whether  the  physics 

embodied  in  the  equations  solved  is  adequate  to  explain  the  simulated 
system. 
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Initially  only  two-dimensional  simulations  of  the  magnetosphere  were 
performed;  LeBoeuf  et  al.  (1978),  Lyon  et  al.  (1980,  1981a).  LeBoeuf  et 
al.  (1978)  used  a  fluid  particle  approach  in  the  MHD  modeling  and  were 
able  to  reproduce  a  Dungey  (1961)  type  magnetosphere  using  the  very 
resistive  Lax  algorithm.  Lax  (1954).  Lyon  et  al.  (1980),  solving  the 
standard  UHD  equations  using  a  combination  of  flux  corrected  transport 
(Boris  and  Book,  1973,  Zalesak,  (1979),  and  partial  donor  cell  method, 
(PDM)  Hain  (1978),  addressed  the  very  dynamic  problem  of,  shock 
focusing.  In  this  work  they  studied  the  effect  of  solar  shocks  on  the 
earth's  magnetosphere.  It  was  found  that  solar  shocks  interacted  with 
the  magnetosphere  in  such  a  way  as  to  focus  the  shock  in  the  tail  region 
at  -15  to  -20  Re«  This  focusing  produced  increased  temperature  and 
pressures  in  this  region.  It  was  speculated  in  this  paper,  as  well  as  in 
Brecht  et  al.,  (1980),  that  such  effects  might  well  trigger  magnetic 
substorms.  In  addition,  Lyon  et  al.  (1981a)  studied  the  problem  of 
substorm-like  events  triggered  by  southward  directed  IMF.  For  this 
problem  a  code  was  developed  that  used  a  leap-frog  time-integration 
scheme  and  20t*1  order  finite-difference  approximation  to  the  spatial 
derivatives  along  with  a  f lux-convected  transport  algorithm  (Boris  and 
Book,  1973,  Zalesak,  1975)  which  provides  the  minimum  possible  numerical 
dissipation  consistant  with  a  solution  showing  no  extraneous  extrema. 
The  results  of  these  simulations  produced  behavior  in  the  tail  which 
tended  to  support  the  neutral-line  model  for  magnetic  substorms  (Dungey, 
1961;  Uones,  1973,  and  McPherron  et  al,  1973).  Research  is  continuing 
into  the  2-D  modeling,  particularly  with  regard  to  resistive  reconnection 
and  how  it  might  be  affected  by  variations  in  the  solar  wind. 
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However,  the  earth's  magnetosphere  is  a  three-dimensional  entity* 
The  need  to  consider  it  as  such  has  been  discussed  by  many  authors,  most 
recently  by  Akasofu  et.  al.  (1978),  Lui  et  al.,  (1978),  Fairfield  (1979), 
and  Coroniti  et  al*  (1980).  These  authors,  to  mention  a  few,  have 
presented  considerable  evidence  that  two-dimensional  models  of  the 
earth's  magnetosphere  are  incomplete.  This  is  in  part  due  to  the 
variability  and  three-dimensional  nature  of  the  solar  wind,  particularly 
the  Interplanetary  Magnetic  Field  (IMF)  carried  by  the  wind.  Within  the 
last  year,  three-dimensional  codes  have  been  constructed  and  the  results 
of  their  initial  work  reported;  (Leboeuf  et  al.,  1981;  Brecht  et  al. , 
1981,  and  Wu  et  al.,  1981).  These  three  simulations  considered  three 
somewhat  different  problems.  The  first  calculation  was  by  LeBoeuf  et  al. 
(1981).  In  this  simulation  a  southward  IMF  interacted  with  a  current 
loop.  The  calculation  used  periodic  bounding  conditions  and  a  particle 
fluid  approach  with  a  Lax  scheme  to  solve  the  MHD  equations.  The 
calculation  was  reported  to  reproduce  a  IXingey  (1961)  like 
magnetosphere.  Unfortunately,  the  scheme  was  so  diffusive  that  close 
inspection  of  the  published  magnetic  field  vector  plots  reveal  no  obvious 
northward  field  such  as  one  expects  to  see  near  the  earth's  dipole 
field.  The  code  also  had  limited  spatial  extent  in  the  y,  z  direction 
causing  the  possibility  of  boundary  effects  to  be  of  concern.  The  strong 
numerical  diffusion  in  addition  to  causing  eradication  of  some,  if  not 
all,  of  the  dipole  field  also  prevents  the  studying  of  any  dynamics  of 
reconnection  regions  in  the  tail. 

Brecht  et  al.  (1981)  studied  the  effects  of  an  east-west  component 
of  the  IMF  being  added  to  a  solar  wind  carrying  a  -2y  southward  field* 
This  calculation  was  fully  three-dimensional  with  no  planes  of  symmetry 


8 


I 


used.  This  was  necessary  to  allow  for  the  rotation  of  the  IMF.  The  code 
used  in-flow  boundary  conditions  on  the  +  x  boundary  and  out-flow 
conditions  on  all  other  boundaries.  The  numerical  scheme  employed  was  a 
partial  donor  cell  method  (PDM),  Hain  (1978).  PDM  couples  a  low  order 
numerical  scheme  to  a  high  order  scheme  in  a  nonlinear  fashion  in  order 
to  prevent  numerical  oscillations  while  maintaining  good  accuracy  and 
minimum  numerical  diffusion.  The  results  showed  that  a  rotation  of  the 
IMF  toward  an  east-west  direction  rotated  or  skewed  the  entire 
magnetotail  and  produced  an  induced  By  in  the  plasma  sheet  of  the  size 
seen  in  the  solar  wind.  These  results  were  found  to  be  consistent  with 
experimental  data;  Hardy  et  al..  (1979),  Lui  (1979),  and  in  agreement  with 
an  analytic  calculation  by  Cowley  (1981). 

The  simulation  by  Wu  et  al.  (1981)  addressed  the  case  where  no  IMF 
exists  in  the  solaT  wind.  They  used  a  Rusanov  numerical  scheme  to  solve 
the  MHD  equations.  The  calculation  was  performed  on  a  quadrant  of  the 
magnetosphere  meaning  that  both  the  meridional  and  equatorial  planes  were 
planes  of  symmetry.  Because  of  the  symmetry  planes  no  reconnection  can 
occur  due  to  numerical  resistivity.  The  boundary  conditions  required 
out-  flow  at  45°  from  the  x  axis.  The  calculation  reproduced,  apart  from 
the  diffusion  effects  noted  by  the  authors,  many  of  the  major  features  of 
the  magnetosphere. 

In  this  paper  we  present  results  of  refined  simulations  where  we 
consider  the  structure  and  dynamics  of  the  magnetosphere  in  the  presence 
of  a  -  2y  IMF.  The  paper  will  discuss  in  detail  the  calculational  set  up 
and  equations  used  in  the  simulation.  In  addition,  we  present  time 
resolved  results  of  island  formation  in  the  tail  and  attendant  flows  and 
densities.  Finally  the  conclusions  derived  from  these  simulations  will 


be  presented 
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Calculation 


As  mentioned  earlier  in  this  paper,  the  code  uses  a  PDM  scheme  for 
solving  the  MHD  equations.  This  scheme  allows  resolution  of  shocks  in  3 
cells  and  therefore  permits  the  use  of. fewer  grid  cells  than  is  necessary 
for  schemes  where  more  numerical  diffusion  is  present,  for  example.  Lax 
or  Rusanov.  The  MHD  equations  the  code  solves  are 


3p 

3t 


V*pv 


«  -  vVv  -  -  (V(P  +  Q)  +  J/c  x  B] 
3t  P 

ap  _  —  2 

~  -  -  V-Pv  -  (y  -  1)  (P  +  Q)V»v  -  nJ 


_  \  _  _  gg 

where  E  -  -  (v  x  B)/c  -  nJ  and  J  *  c7x  B/4*  -  —  . 

The  variables  p,v,P,B,n,  and  Q  are  mass  density,  velocity,  p  asraa 
pressure,  magnetic  field,  resistivity  and  artificial  viscosity. 
Artificial  viscosity  is  necessary  to  insure  that  the  Rankine-Hugoniot 
relations  are  satified  throughout  the  code.  It  is  only  activated  in  the 
region  of  the  simulation  where  there  is  compressive  flow,  7»v  <  0.  The 
PDM  algorithm  in  conjunction  with  this  set  of  equations  guarantees  that 
numerically  V»B  ■  0  everywhere  in  the  code,  a  critical  feature  to  the 
validity  of  the  simulation  not  automatically  obtained  by  all  schemes. 
The  effect  of  high  Alfven  speeds  in  requiring  unreasonably  short  time 
steps  is  controlled  by  the  3¥/3t  term  in  Amperes'  law.  The  inclusion  of 
the  displacement  current  effectively  limits  the  propagation  speed  of  the 
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Alfven  waves  Co  Che  speed  of  lighc,  which  can  be  arbiCrarily  seC  in  Che 
code*  The  boundary  condiCions  for  Chis  simulation  are  Che  same  as  those 
used  in  our  previous  3-D  paper  [BrechC  eC  al.,  1981]*  The  +  X  boundary 
has  an  inflow  condiCion  and  all  other  boundaries  have  an  ouCflow 
condition. 

The  simulation  presented  in  this  paper  was  performed  on  a  grid  using 

a  stretched  mesh  of  29  x  21  x  21  cells  in  the  x,y  and  z  directions  of  the 

solar-magnetospheric  coordinate  system.  The  spatial  extent  of  the 

numerical  grid  is  -  100  R  <  X  <  30  R  ,  -  40  R  <  Y  <  40  R  and  -  40  L  < 

e  e  e  e  e 

Z  <  40  Re.  The  earth's  dipole  field  is  formed  by  placing  a  current  loop 
atX»V*Z»0  with  the  appropriate  current  to  replicate  the  earth's 
field.  A  solar  wind  is  then  introduced  into  the  calculation  from  the  +  x 
boundary,  which  distorts  the  dipole  field  and  creates  the 
magnetosphere.  The  solar  wind  may  carry  an  IMF  of  any  strength  and 

orientation.  The  solar  wind  introduced  into  the  calculation  presented 

.j  —3 

here  has  a  velocity  of  400  km  sec  ,  a  density  of  5  cm  and  a  Bz  * 

_  O 

-  2y.  The  grid  was  initially  loaded  with  a  plasma  density  of  1  cm  . 
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Results 


The  results  presented  here  are  a  3-D  time  evolution  calculation  of 
the  magnetosphere  as  a  southward  (-2f)  IMF  Interacts  with  the  earth's 
magnetic  field.  Reconnection  was  permitted  to  occur  In  this 

simulation.  We  had  the  option  of  no  reconnection,  using  real  resistivity 
in  the  equation,  or  permitting  numerical  reconnection.  We  chose  the 
latter  case  for  several  reasons.  Foremost  among  these  is  the  fact  we  did 
not  know  where  to  put  the  resistivity  and  in  what  amount.  Arbitrary 
prescriptions  did  not  appeal  to  us.  Numerical  resistivity  occurs  mainly 
where  opposed  fields  are  being  driven  together  in  a  compressive  fashion 
causing  the  location  of  reconnection  regions  to  be  determined  by  the 
dynamics  of  the  calculation  not  the  predisposition  of  the  author.  At 
this  juncture  of  our  research  we  are  actively  testing  various  ways  to 
model  the  microinstabilities,  Lyon  et  al.  (1981b),  which  are  thought  to 
provide  the  necessary  resistivity  for  the  reconnection  process.  However 
the  present,  the  proper  physics  and  its  implementation  is  not  well 

resolved.  It  should  be  noted  that  for  steady  state  configurations  such 
as  those  used  in  the  research  of  Birn  and  Hones  (1981)  the  PDM  scheme 
will  not  numerically  reconnect.  The  numerical  reconnection  shown  here  is 
purely  a  driven  phenomenon.  One  disadvantage  of  using  numerical 

reconnection  is  that  with  this  code  there  is  no  resistive  heating  of  the 
plasma.  Therefore,  only  the  magnetic  field  and  pressure  forces  caused  by 
the  compressive  flow  will  affect  the  plasma  dynamics.  As  will  be  shown 
in  this  paper  these  terms  alone  provide  rather  Impressive  effects  in  the 
earth's  magnetotail. 

The  calculation  begins  with  the  earth's  dipole  field  being  acted 
upon  by  a  solar  wind  of  the  previously  stated  properties  as  well  as 
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the  -  2y  southward  IMF.  A  reconnection  point  forms  in  the  near  earth 
tail  region  approximately  4500  secs  into  the  calculation.  (In  this  paper 
all  times  quoted  are  real  time  and  not  computer  simulation  time).  This 
reconnection  point  begins  to  move-  tailward.  Eventually  another 
reconnection  point  occurs  and  a  magnetic  island  forms  and  begins  to  move 
in  tailward  direction.  As  this  island  moves  tailward  another  island  is 
formed.  The  general  dynamics  of  the  system  is  quite  similar  to  the  2-D 
simulations  reported  by  Lyon  et  al.  (1981a)  and  the  3-D  simulation  by 
Birn  and  Hones  (1981).  One  main  difference  in  the  results  of  the  2-D  and 
3-D  simulations  is  that  for  the  3-D  simulations  strong  earthward  flows 
are  observed  whereas  in  2-D  this  was  not  the  case.  It  was  speculated  by 
Lyon  et  al.  (1981a)  that  the  lack  of  earthward  flow  could  be  due  to  the 
two-dimensionality  of  the  calculation  and  this  appears  to  be  the  case. 

In  figures  1  -5  the  global  temporal  dynamics  of  the  various  plasma 
parameters  are  shown.  Each  set  of  figures  begins  at  7691  sec  (0  min) 
into  the  calculation  and  run  through  11351  secs  (~  65  min)  and  therefore 
cover  slightly  over  one  hour  of  real  time.  The  beginning  time  was 
selected  because  most  of  the  initial  transients  produced  as  the  dipole 
field  was  transformed  into  a  magnetospheric  topology  were  out  of  the 
system.  One  remnant  of  such  a  transient  is  seen  as  a  reversed  field 
structure  in  the  tall  of  the  magnetosphere  in  Fig.  1  at  about  55  Re.  In 
addition  we  have  chosen  two  locations  in  the  magnetotail  at  Y  ■  0,  Z  »  0, 

and  X  -  19  Re  and  -  32  Rg  to  display  plasma  properties  that  a  satellite 

would  detect  during  this  time  frame  (Figs.  6  -  10). 

Fig.  1  is  a  time  sequence  of  magnetic  vector  plots  on  the  x-z 
meridional  plane.  The  vectors  are  normalized  to  the  local  magnitude  of 
B,  so  the  size  of  the  vector  provides  no  indication  as  to  the  strength  of 
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Fig.  1  —  Time  evolution  of  the  Bx  and  Bs  magnetic  field  component*  plotted  in  the 
meridional  (noon-midnight)  plane.  For  purposes  of  discussion  7698  and  11351  sec 
correspond  to  0  min.  and  ~  65  min.  respectively.  The  length  of  the  vector  drawn 
provides  no  information  about  the  magnitude  of  the  fields.  All  distances  are  mea¬ 
sured  in  earth  radii. 


the  magnetic  field.  In  Fig.  la  one  notes  the  existence  of  a  reconnection 
point  at  about  -68  R0  in  the  tail  as  well  as  an  o-point  in  the  -35  Rg  to 
-  40  Rg  region.  In  the  following  several  snapshots  one  sees  the  interior 
island  grow  in  the  tailward  direction  and  by  9766  secs  is  beginning  to 
accelerate  tailward.  The  last  timestep  shows  the  formation  of  a  new 
island  with  the  inner  x-point  located  at  approximately  -20  Re  in  tail. 

The  velocities  that  correspond  to  the  creation  and  motion  of  these 
sets  of  neutral  lines  or  o-points  can  be  seen  in  Fig.  2.  Here  we  show 
the  velocities  in  the  equatorial,  x-y  plane.  In  these  plots  the  arrow 
length  indicates  the  magnitude  of  velocity.  The  unit  vector  is  set  to  be 
400  km  sec  the  solar  wind  velocity.  The  general  features  of  the 
magnetosphere  can  be  see  in  Fig.  2a.  The  stagnation  point  in  the  bow  is 
located  at  about  X  =  12  Re,  and  the  magnetopause  can  be  seen  draped 
around  to  the  sides.  In  the  interior  region  of  the  magnetotail  one  notes 

a  demarcation  line  of  very  low  velocities  at  X  -  29  Rg  corresponding  to 

the  location  of  the  x-point  seen  in  the  magnetic  field  plots  at  the  same 
time.  Earthward  of  this  point  there  is  flow  toward  the  earth.  Tailward 
of  this  point  the  flow  is  predominately  away  iron  the  earth.  It  should 
be  noted  that  in  the  meridional  plane  one  sees  these  flow  patterns 
restricted  primarily  to  -  2  Se  <  Z  <  2  R  g.  To  a  large  extent  the 
magnitude  of  the  flow  is  determined  by  the  mass  loading  of  the  magnetic 
field.  The  correspondence  between  high  flow  regions  and  low  densities  is 
quite  high  as  can  be  seen  in  Fig.  2  and  4,  which  show  the  perspective 
plots  of  velocity  and  density  in  the  equatorial  plane. 

As  one  scans  through  the  various  time  steps  shown  in  Fig.  2,  one 
notes  the  velocities  pulsate  within  a  period  of  about  20-30  minutes. 
This  is  noticeable  if  one  considers  the  X  ~  -  22  Rg  grid  line  shown  in 
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Fig.  2  —  Time  evolution  of  Vx  and  Vy  velocities  in  the  equatorial  plane.  The  unit 
vector  is  400  km  sec-1  which  is  the  sohur  wind  velocity. 


Fig.  2.  As  time  progresses  one  notes  two  obvious  features.  The  first  is 
that  the  speed  of  the  flows  increases  with  time  as  successive  islands  are 
formed.  This  feature  was  noted  in  the  2-D  simulation  by  Lyon  et  al. 
(1981a).  In  both  calculations  each  succeeding  reconnection  event  reduces 
the  amount  of  mass  on  the  field  lines  and  consequently  the  flow  speed 
becomes  larger.  Fig.  2f  also  shows  the  second  feature,  an  asymmetry  in 
the  flew  pattern.  This  is  a  particularly  interesting  feature  caused  by 
the  motion  of  the  x-point  at  the  bow  of  the  magnetosphere.  Because  of 
the  large  number  of  timesteps  in  this  calculation  numerical  noise 

permitted  the  reconnection  point  on  the  bow  to  shift  slightly.  Once  this 

occurs  the  continued  flow  begins  to  develop  an  asymmetry  on  the  bow  which 

manifests  itself  as  motion  in  the  tail  and  produced  the  asymmetic  flow, 
much  as  one  might  expect  in  the  real  magnetosphere.  The  density  and 
pressure  are  noted  to  be  asymmetric.  In  this  simulation  the  asymmetry 
existed  in  only  the  dawn  -  dusk  direction.  The  fact  that  the 
reconnection  point  has  become  slightly  asymmetric  may  be  more  realistic 
than  the  usual  perfectly  symmetric  simulation.  It  should  also  be  noted 
that  the  region  of  very  rapid  flow  is  restricted  mainly  to  the  equatorial 
plane,  but  as  can  be  seen  does,  in  fact,  cover  a  large  spatial  region  in 
both  the  x  and  y  directions.  The  flow  speeds  reach  velocities  as  high  as 
650  km  sec-1  in  the  earthward  direction  and  800  km  sec  *  tailward. 

Coincident  with  the  cycling  of  the  island  formation  one  can  see 
considerable  effects  on  the  density  and  pressure  in  the  magnetotail. 
Looking  at  the  density  in  the  meridional  plane  (Fig.  3)  one  sees  two 

important  features.  The  first  is  that,  as  the  successive  neutral  lines 
move  tailward  the  plasma  sheet  thins  tremendously,  and  second  that  the 
plasma  sheet  density  is  considerably  reduced.  In  Fig.  4f  the  density  in 


12 


Fig.  3  —  Perspective  plots  of  density  in  the  meridional  plane.  Units  for  the  density  scale 
are  in  particle  cm-3.  The  rather  jagged  appearance  of  the  bow  shock  is  due  to  the  per¬ 
spective  plotter,  not  data.  This  is  true  for  all  perspective  plots  shown. 


Fig.  4  —  Perspective  plots  of  density  in  the  equatorial  plane 


the  middle  of  the  tail  is  quite  low,  but  an  enhancement  of  the  density  is 
occurring  on  the  flanks  of  the  tail.  The  asymmmetry  caused  by  the 
shifted  reconnection  point  is  only  about  10%  in  density  from  the  dawn 
side  to  the  dusk  side.  Fig.  5  shows  the  pressure  in  the  meridional  plane 
where  one  sees  the  rather  dramatic  sheet  thinning  often  associated  with 
reconnection  events,  Hones  and  Schindler  (1979),  and  Hones  (1979).  In 
Fig.  5a-f  one  notes  the  progressive  thinning  of  the  sheet  as  indicated  by 
the  pressure  dropout  and  a  general  overall  reduction  in  pressure  along 
the  plasma  sheet,  Z  =  0.  Also,  one  notes  a  slight  recovery  of  the 
pressure  occurring  inside  of  about  X  =  -  20  Rg  and  moving  tailward. 

In  the  next  several  paragraphs  we  discuss  data  taken  from  the  code 
at  two  locations  in  the  magnetotail.  The  time  resolution  of  the  data  is 
approximately  5  min  intervals.  In  a  sense  this  is  plasma  and  field  data 

as  a  satellite  would  detect  it  during  the  hour  we  have  displayed.  The 

two  locations  are  located  at  Y  ■  0,  Z  =  0,  X  ~  -  19  Re  and  -  32.  Rg. 

In  Fig.  6  we  show  the  x  velocity  as  seen  at  the  two  locations 
mentioned  above.  The  inner  location  sees  an  earthward  flow  which  becomes 
progressively  stronger  with  time.  There  is  a  slight  oscillation  in  the 
velocity,  but  generally  it  is  increasing.  The  maximum  velocity  is 
reached  at  about  55  min.  At  this  time  the  magnetotail  neutral  line  has 
moved  to  beyond  the  X  ~  -32  Rg  location.  Fig.  le.  Shortly  after  this 

time  the  tail  reconnects  and  a  x-point  forms  at  about  X  ~  -22  Re,  as  can 

be  seen  on  the  Bz  field  evolution  in  Fig.  If.  The  creation  of  the  inward 
neutral  line  corresponds  to  the  rather  precipitous  drop  in  velocity  seen 
at  the  inner  point.  The  outer  location  experiences  an  oscillating 
velocity  of  about  100  km  sec”*.  During  the  outward  stretching  phase  at 
about  55  min  the  tailward  flow  drops  to  about  50  km  sec”*.  Locations 
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Fig.  5  —  Perspective  plots  of  pressure  in  the  meridional  plane, 
Units  for  the  pressure  scale  are  ergs  cm-3. 


slightly  closer  to  the  earth  experience  a  complete  flow  cycle  from 
negative  to  positive  and  back  to  negative  flow.  As  soon  as  the  x-point 
forms  earthward  of  the  observation  point  one  sees  the  tailward  velocity 
go  from  -50  km  sec  *  to  -650  km  sec  ^  in  less  than  5  min.  Similar  to  the 
results  of  Lui  (1979)  and  Hones  and  Schindler  (1979),  strong  dawn-dusk 
velocities  are  detected  at  these  two  locations.  The  inner  point 
displayed  primarily  a  duskward  flow  until  the  neutral  line  appeared  at 
-22  Re«  At  this  time  the  flow  changed  rapidly  to  -300  km  sec-*  in  the 
dawnward  direction.  The  outer  point  in  the  meantime  displayed  a 
constantly  increasing  duskward  flow  which  peaked  at  -160  km  sec-*.  The 
location  of  the  observation  point  with  respect  to  the  neutral  line  seems 
to  determine  which  direction  the  dawn-dusk  flows  go,  Fig.  7.  The  z 
velocity  generally  stayed  under  10  km  sec-*  for  both  locations  and 
oscillated  between  the  positive  and  negative  direction. 

The  Bx  field  was  quite  variable  with  the  inner  point  experiencing 
larger  fields  as  expected.  The  B^  field  in  these  simulations  remained 
negligible  thoughout  the  simulation  as  detected  at  these  locations.  In 
Fig.  8,  the  time  evolution  plots  of  the  z  magnetic  fields  are 
displayed.  Here  one  sees  some  very  interesting  behavior.  The  inner 
point  experiences  fairly  large  positive  Bz  components  ranging  from  5  to  8 
gamma.  However,  at  the  time  the  neutral  line  is  estabilished  at  ~  -  22 
Re  one  sees  a  precipitous  drop  in  the  Bz  component.  This  is  caused  by 
the  close  proximity  of  the  neutral  line.  The  outer  point  experiences 
more  dynamic  changes  in  Bg.  Initially  the  neutral  line  is  earthward  of 
this  point.  One  can  correlate  the  expansion  of  the  magnetic  island  seen 

in  Fig.  la-c  with  the  general  decrease  of  Bz  from - ly  to  ~  -  1.25y. 

The  Island  then  moves  tailward  and  passes  beyond  the  outer  point  as  shown 
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Fig.  7  —  Temporal  variations  of  Vy  as  seen  by  the  two  tail  ward  locations. 
Squares  indicate  X  =  -  19  Re  and  triangles  X  =  -  32  Re. 
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Fig.  8  —  Temporal  variations  of  Bz  as  seen  by  the  two  tail  ward  locations. 
Squares  indicate  X  *  -  19  Re  and  triangles  X  *  -  32  Re. 
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by  the  change  of  Bz  from  ~  -  1.25y  to  ~  1.5y.  As  discussed  earlier  and 
coincident  with  the  rapid  increase  of  plasma  flow  a  neutral  line  is 
established  earthward  of  this  point  and  the  Bz  changes  from 
~  1.5y  to  - 2y  in  about  the  same  5  minute  interval. 

The  relationship  between  the  pressure,  density  and  Bz  component  of 
the  magnetic  field  is  somewhat  confused.  Figure  9  shows  that  the  two 
locations  in  the  tail  experience  approximately  the  same  pressure 
variation  with  perhaps  a  slight  time  delay.  However,  the  formation  of 
the  new  x-point  causes  the  pressure  behavior  to  change  at  about  65  min. 
At  this  time  both  locations  show  quite  similar  pressure  levels.  The 
inner  location  shows  a  drop  of  a  factor  of  approximately  two  and  the 
outer  point  is  enhanced  by  about  the  same  factor.  In  the  density  plot. 
Fig.  10,  a  build  up  of  density  occurs  at  about  20  min.  Subsequent  to 
this  a  general  decline  in  density  occurs.  Hie  build  up  is  associated 
with  the  recovery  of  the  plasma  sheet  from  the  passage  of  the  previous 
island  now  located  far  in  the  tail,  Figs,  lc  and  d.  The  decline  in 
density  after  the  20  min  mark  is  associated  with  the  movement  tailward  of 
the  o-point  located  at  X  ~  -  35  Rg,  see  Fig.  Id  and  e,  and  the 

associated  stretching  out  of  the  tail.  However,  as  can  be  seen  the 

density  plots  do  not  display  the  variations  seen  in  the  pressure  plots. 

Comparing  the  pressure  plots,  Fig.  9,  with  magnetic  field  plots, 

Fig.  8,  one  notes  that  the  outer  location  shows  increases  in  pressure 
associated  with  decreases  in  the  magnitude  of  Bz.  This  relationship 
holds  until  about  the  60  min  mark,  by  this  time  an  x-point  has  passed  by 

this  location  and  a  new  one  is  about  to  form  at  X 22  Rg.  The 

pressure  is  found  to  increase  by  a  factor  of  about  two.  The  magnitude  of 
the  Bz  field  passes  through  zero  and  ends  up  larger  than  it  has  been 
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Fig.  9  —  Pressure  as  a  function  of  time  for  the  location 
X  *  -  19  Re  (squares)  and  X  *  -  32  Re  (triangles) 
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Fig.  10  —  Density  as  a  function  of  time  for  the  locations 
X  =  -  19  Re  (squares)  and  X  =  -  32  Re  (triangles) 


23 


during  this  whole  data  sequence .  The  inner  point  exhibits  different 
behavior.  For  the  first  20  min  the  pressure  and  the  magnitude  of  Bz 
follow  one  another.  Then  the  pressure  drops  while  Bz  increases.  During 
the  final  30  min  of  data  the  two  quantities  then  follow  one  another 
again.  The  correspondence  of  the  density  to  either  pressure  or  Bz  is 
weak.  Only  the  build  up  of  density  at  20  min  appears  to  effect  the 
variations  seen  in  either  the  pressure  or  Bz» 

These,  then,  are  the  basic  results  from  the  calculation.  In  the 
next  section  we  discuss  implications  of  these  results  and  how  they 
correlate  with  experimentally  obtained  data. 
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Discussion 


While  one  encounters  considerable  difficulty  in  correlating  in  an 
exact  fashion  the  published  experimental  results  with  those  of  a 
numerical  simulation,  the  experimental  results  do  provide  a  guideline  to 
the  kind  of  questions  one  ought  to  ask  concerning  the  numerical 
results.  There  have  been  many  articles  published  reporting  satellite 
data  from  all  areas  of  the  magnetosphere  during  substorm  events. 
However,  for  the  purposes  of  this  discussion  we  have  restricted  ourselves 
to  the  data  reported  fr^tn  the  IMP-6  satellite  by  Hones  and  Schindler 
(1979),  and  Lui  (1979).  It  is  not  our  intention,  nor  within  our 
capability,  to  resolve  differences  between  their  respective  analyses  of 
the  IMP-6  data.  Our  intention  here  is  to  compare  the  reported  data  with 
the  calculation  we  have  performed. 

As  discussed  earlier  we  have  selected  and  shown  temporal  data  from 
various  locations  in  the  magnetosphere.  The  detailed  data  sho'-bn  by  Hoi..*^ 
and  Schindler  (1979)  are  taken  primarily  in  the  X  ~  -  30  Rg  region,  but 
varies  somewhat  in  the  y  and  z  direction.  Lui  (1979)  displays  detailed 
data  in  the  same  region,  but  shows  not  only  data  at  X  ~  -  30  Re,  but 
X  ~  -  20  Rg  as  well.  In  their  statistical  analysis  the  respective 
authors  included  data  from  all  positions  in  x,  subject  to  the  various 
criteria  placed  on  the  data  to  determine  its  applicability. 

At  the  risk  of  over-  simplification  of  these  authors'  results  one 
can  sum  up  the  major  results  in  the  following  way.  As  a  general  rule  the 
measured  flows  appeared  to  be  field  aligned.  There  appeared  coincident 
with  these  events  a  dawn-dusk  flow  favoring  the  duskward  direction.  Lui 
(1979),  arguing  in  favor  of  the  current  interruption  model  which  we 
cannot  assess  with  this  particular  calculation,  reported  that 
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statistically  speaking  there  was  not  a  strong  correlation  between 
negative  Bz  and  tailward  flow.  Hones  and  Schindler  (1979)  found  that  for 
"40  -  50  percent  of  the  events,  no  substorm  effect  on  the  flow  is  seen: 
in  another  40  -  50  percent  there  is  tailward  flow  associated  with  the 
substorms;  and  in  10  -  15  percent  there  is  earthward  flow".  In  addition, 
to  comparing  AE  index  to  plasma  flows,  both  Lui  and  Hones  and  Schindler 
also  correlated  plasma  flow  and  sheet  thinning.  The  respective  authors 
reached  apparently  opposite  conclusions  based  on  their  IMP-6  data  for  the 
relation  between  plasma  sheet  thinning  and  flow  direction  during  a 
substorm. 

Our  calculations  do  not  resolve  the  differences  of  these  authors. 
For  one  thing  it  is  not  clear  that  the  solar  wind  conditions  in  our 
simulation  accurately  represent  the  solar  wind  conditions  that  were 
present  during  the  events  for  which  t'^se  authors  show  data.  As  we  have 
found  earlier,  [Brecht,  et  al.  (1981)],  inclusion  of  effects  such  as  By 
in  the  solar  wind  makes  considerable  modifications  to  the  magnetotail, 
not  to  mention  the  possible  effects  on  microinstabilities  thought  to 
provide  resistivity  for  tearing.  We  can,  however,  from  our  simulation 
either  add  to  the  confusion  or  provide  some  clarification,  depending  on 
the  reader's  point  of  view. 

It  is  heartening  to  note  that  the  magnitude  of  the  quantities 
displayed  in  Figs.  6-10  are  similar  to  those  found  in  the  magnetotail. 
We  see  in  Fig.  8  variations  and  magnitudes  of  Bz  consistent  with  the  data 
shown  by  Lui  (1979),  see  Figs.  14  and  16,  which  correspond  with  our 
selected  observation  data  points.  The  velocities  seen  in  our  calculation 
are  also  consistent  with  those  seen  during  substorm  events.  The 
pressures  also  appear  to  be  consistent  with  the  data  albeit  somewhat 
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lower,  possibly  due  Co  Che  lack  of  ohmic  hearing  expected  Co  occur 
because  of  resistivity.  We  do  noC,  however,  see  dropouts  of  several 
orders  of  magnitude  in  the  data  at  Y  »  0,  Z  *  0.  Factors  of  five  seem  to 
be  the  most  we  see  directly  in  the  magnetotail.  However,  consideration 
of  Figs.  3-5  show  that  any  displacement  from  the  central  part  of  the 
plasma  sheet  either  in  y  or  z  can  result  in  considerable  drop  in  either 
pressure  or  density. 

We  have  found  in  this  simulation  that  the  data  one  collects  is 
extremely  dependent  on  location.  Consider  Fig.  2  at  a  position  of  Z  *  0, 
Y  ~  >  6  Ke  and  X  between  -  23  R£  and  -  45  Rg.  One  notes  that  if  the 
satellite  had  been  located  here  during  the  hour  that  several  x-points 
existed  in  the  tail  virtually  no  flow  modification  would  have  been 
detected.  This  type  of  event  was  reported  by  Hones  and  Schindler  (1979) 
at  Y  *  10.6  Re.  However,  at  Y  <  -  6  Re  during  this  sane  tine  one  would 
have  seen  flows  approaching  400  km  sec-^.  At  the  point  X  =  -  22.8  Re  the 
displacement  in  y  also  changes  the  characteristics  of  the  flow  because  of 
curvature  of  the  neutral  line  itself.  An  additional  feature  of  this 
calculation  is  the  dependence  on  Z  seen  in  terms  of  earthward  and 
tailward  flows.  All  the  flows  shown  here  are  on  the  equatorial  plane,  - 
2  Rg  <  Z  <  2  Re.  If  one  moves  above  or  below  the  Z  *  0  plane  one  quite 
frequently  finds  flows  of  smaller  magnitude.  These  flows  occurred  in  x 
locations  where  strong  flows  existed  at  Z  «  0. 

Simultaneous  data  taken  at  the  two  locations  in  the  sheet  show  the 

effects  of  the  x-point  motion.  In  Fig.  8  one  sees  that  the  point  at 

1 

X  -  19  Kg  experiences  variations  in  Bz,  which  remains  positive  although 

it  drops  in  magnitude  as  the  x-point  forms  close  to  it.  The  outer  point 
sees  a  reversal  in  Bz  as  the  x-point  passes  this  location  and  then 
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returns  to  negative  as  the  new  x-point  forms  at  ~  22  Rfi.  Looking  at  Fig. 
6  one  sees  that  the  inner  point  detects  earthward  flow  most  of  the  time, 
but  the  outer  point  sees  only  tailward  flow,  although  the  sign  of  Bz 
reversed.  As  the  new  x-point  forms  the  tailward  flow  jumps  to  over  -  600 
km  sec-*.  Not  shown  in  these  figures  is  a  point  located  at 

X  ~  -  27  Re«  Here  an  observer  would  not  only  experience  reversal  of  Bz, 

but  also  reversal  of  the  flow  pattern  from  tailward  to  100  kra  sec  * 
earthward  back  to  >  700  km  sec-*  tailward.  The  point  to  be  made  is  that 
during  the  passage  of  the  x-point  the  dynamics  of  the  new  x-point 
formation  caused  the  data  to  be  extraordinarily  sensitive  to  location. 
The  existence  of  dusk-dawn  flows  were  also  found  at  our  two  data 
points.  The  flow  was  duskward  for  most  of  the  simulation,  but  when  the 
new  x-point  formed  the  outer  location  experienced  a  tremendous  change. 

The  flow  went  from  140  km  sec-*  to  -  350  km  sec-*.  The  pressure  shows 

the  type  of  plasma  sheet  dropouts  seen  in  the  data.  However,  one 
interesting  point  should  be  noted.  As  the  neutral  line  approached  the 
inner  point  at  65  min. ,  the  plasma  pressure  dropped  along  with  the  plasma 
flow  speed  and  the  magnitude  of  Bz«  At  the  outer  point  the  opposite 
result  occurred  at  this  time. 

Finally,  as  reported  by  Lui  (1979),  and  Hones  and  Schindler  (1979) 
we  find  that  most  of  the  energetic  flows  appear  to  be  field  aligned. 
These  flows  are  found  to  originate  at  the  neutral  line.  In  fact,  they 
cannot  be  generated  eleswhere  because  n  ■  0  and  only  J  x  B  and  pressure 
forces  are  in  effect.  A  comparison  of  Figs.  1  and  2  illustrates  this 


point.  While  Fig.  1  does  not  display  magnetic  field  lines,  By  is 
typically  an  order  of  magnitude  smaller  than  the  other  components. 


In  short,  within  this  one  simulation  we  find  conditions  and  plasma 


behavior  quite  similar  to  the  various  statistical  features  reported  by 
Hones  and  Schindler  (1979),  and  lui  (1979).  Whether  or  not  they  would  be 
detected  depended  entirely  on  the  observer's  locations  in  the  magnetotail 
at  the  time  of  the  substorm  like  event* 


Conclusions  and  Summary: 

The  3-D  time  dependent  simulations  of  the  magnetosphere  presented  in 

this  paper  display  many  of  the  features  commonly  associated  with 

substorras.  It  was  found  that  motion  of  the  x-point  caused  plasma  dropout 

to  occur  in  the  magnetotail.  The  flows  seen  in  this  calculation  appeared 

to  be  oriented  in  a  field  aligned  direction.  The  speed  of  flow  was 

inversely  related  to  the  plasma  density  on  the  field  line.  As  the 

simulation  progressed  the  density  declined  in  the  tail  and  the.  flow 

speeds  increased.  The  magnitude  of  these  flows  and  of  other  plasma  and 

magnetic  parameters  appears  consistent  with  satellite  data  during 

substorm  periods,  however,  until  the  mass  on  the  field  lines  is  removed, 

the  reconnection  events  do  not  produce  behavior  that  corresponds  to  the 

usual  picture  of  substorms.  This  fact  is  reasonable  because  the  local 

Alfven  speed  is  a  good  estimate  for  the  maximum  velocities  produced  by 

reconnection.  tor  plasma  parameters  typically  found  in  the 

_3 

tail,  B  ~  2U  y  and  n  ~  1  cm  the  Alfven  speed  is  400  km/sec.  In  the 
neutral  sheet  the  Alfven  speed  drops  to  40  km/sec.  Therefore,  in  order 
to  have  large  scale  accelerations  one  needs  plasma-dropout  to  achieve  the 
500  -  1000  km  sec-*  speeds  reported  for  substorms.  Thus,  if  the  neutral 
line  model  for  substorms  is  correct,  our  results  indicate  that 
reconnection  Is  a  necessary,  but  not  sufficient  condition  for  jubstorm 
formation. 

The  simulation  exhibited  the  tendency  for  multiple  x-points  or 
neutral  lines  in  the  tail.  This  result  is  consistent  with  Lyon  et.  al., 
(1981).  The  inner  x-point  for  the  time  segment  shown  here  tended  to 
remain  at  X  ~  -  27  Rg  in  the  tail  until  the  o-point  grew  sufficiently. 
It  then  began  to  move  tailward,  stretching  the  tail.  At  this  point  in 
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time  a  new  x-point  was  seen  to  form  at  X  ~  -  20  Rg.  An  interesting 
aspect  of  the  simulation  was  the  temporal  behavior  at  various  locations 
in  the  tail.  Due  to  the  multiple  x-point  formation  and  the  asymmetric 
flow  caused  by  the  x-point  shift  on  the  bow,  it  was  found  that  almost  all 
of  the  various  categories  of  events  reported  by  Lui  (1979),  and  Hones  and 
Schindler  (1979)  could  be  observed  in  this  one  simulation.  The  location 
that  remained  earthward  of  x-points  experienced  no  field  reversal  nor  a 
flow  reversal  in  vx  as  the  various  x-points  were  formed  and  moved 
tailward.  Points  shifted  in  y  from  Y  =  0  often  displayed  no  modification 
to  plasma  flow  during  these  events.  Yet  as  the  flow  moved  toward  the 
dawn  side  due  to  the  shift  in  the  bow  reconnection  point  one  found  that 
at  Y  <-  6.7  Re  flows  were  detected  to  coincide  with  the  formation  of  the 
x-point  while  at  Y  >  6.7  Ke  no  effects  were  detected.  Observation  points 
were  found  in  the  tail  at  Y  *  0,  Z  -  0,  X  ■  -  27  Rg  and  -  32,  Re  which 
showed  field  reversal  from  negative  Bz  to  positive  Bz  and  a  return  to 
negative  as  a  neutral  line  passed  these  locations  and  a  new  one  was 
formed  closer  to  the  earth.  These  points  experienced  different  flow 
behavior  with  the  reversal  of  Bz.  The  point  located  at  X  =  *  27  Rg 
experienced  reversals  in  the  vx  direction  consistent  with  the  neutral 
point  being  either  earthward  or  tailward  of  that  point.  The  point  at  X  = 
-  32  Re  did  not  experience  reversal  of  vx,  but  only  saw  a  reduction  of 
the  tailward  flow.  This  lack  of  the  flow  reversal  can  probably  be 
attributed  to  the  formation  of  the  x-point  at  X  ~  -  20  Re  within  five 
minutes  of  the  previous  neutral  line  passage. 

Finally,  the  various  statistical  categories  used  by  Hones  and 
Schindler  (1979)  and  Lui  (1979)  to  categorize  plasma  and  magnetic 
behavior  during  the  many  substorm  events  experienced  by  IMP-6  appear  to 
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be  found  in  this  one  simulation*  This  leads  one  to  the  not  so  profound 
conclusion  that,  because  of  the  complicated  magnetotail  behavior,  spatial 
and  temporal  variations  allow  one  to  sample  a  full  range  of  dynamic 
behavior  determined  primarily  by  the  location  of  the  observer. 
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